Numerical Extraction of Distributions of Space-charge and Polarization from Laser 

Intensity Modulation Method 
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The Fredholm integral equation of tlie laser intensity modulation method is solved with the 
application of the Monte Carlo technique and a least-squares solver. The numerical procedure is 
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Space-charge and polarization distributions in materi- 
als can be determined by the Laser Intensity Modulation 
MethodfLIMM), which was proposed by Lang and Das- 
Gupta y [see Ref. 0, for a recent review on the LIMM] . 
In the LIMM experiment, a thin sample with electrodes 
on both surfaces is irradiated with a modulated laser. 
The light is absorbed at the irradiated electrode and the 
heat diffuses into the sample. Space-charge and polar- 
ization in the material respond to a non-uniform tem- 
perature distribution, and produce a pyroelectric current 
density j. The real jr and imaginary ji parts of the cur- 
rent are recorded for each modulation frequency v. The 
typical frequency range is between 10 Hz and 100 kHz. 
Lang and Mellinger (21] have recently illustrated how 
to reconstruct the space-charge and polarization profiles 
with polynomial and TikhonovQ regularizations, respec- 
tively. In this letter, an alternative numerical method is 
presented and tested. 

The measured current density_in LIMM is expressed 
as a Fredholm integral equation 



ym 



ji^) = jr + iji = I ItTVS" 



r{z)1{v,z)Az (1) 



where g(z) is the unknown distribution of space charge 
or polarization, ^{v, z) is the solution of one-dimensional 
heat-conduction equation, z is the coordinate in the 
thickness direction, and i — \f—\. The integral in Eq. 
is evaluated over the samples thickness s. For a free 
standing film, T(j^, z) is|3,| 



1{v, z) = A(k/3)-^ cosh[/3(s - z)]csch(/3s) 



(2) 



Here, A is a constant, k is the thermal conductivity 
of the sample. /3 is the complex thermal wave [/3 = 
(1 -|- i)vTruiy^ and D is the thermal diffusivity of the 
material] . 

A numerical method based on the Monte Carlo method 
have been proposed by Tuncer and Gubahski jfl] to solve 
integral equations in the form of Eq. (Q). The method 



has previously been applied to extract the distribution of 
relaxation times from dielectric spectroscopy data [6|, Q 
and the spectral density function (distribution) of dielec- 
tric mixtures In those problems, the distributions 
were probability densities of relaxation times and spec- 
tral parameters. Therefore, a constrained least-squares 
algorithm was implemented which yielded only positive 
valued outputs. In the present problem, the space-charge 
and the polarization distributions g can have both pos- 
itive and negative magnitudes. As a result, we employ 
a similar numerical procedure as before |^. However, 
here a linear least-squares solver is adopted without any 
a-priori assumptions instead of the constrained least- 
squares method. 

The numerical procedure proposed is as follows, 

1. Eq. 1^ is written as a summation over some num- 
ber n of z values for each experimental frequency 
point Vm (w = 1, • ■ ■ , M, M is the number of ex- 
perimental data points). 



j{Vrn) = « 2TrVmS ^ ^ g(z„)T(z^„ 



(3) 



n<M 



2. The Zn values are preselected randomly from a lin- 
ear distribution, < z„ < s. 

3. Eq. Q is written as a matrix equation, 

Tg=j, (4) 

and solved for g. In Eq. Q), T is the mxn 'kernel' 
matrix, where T = T(j^m, Zn), g is a column vector 
with length n (g = g„), and j is a column vector 
with length to. 

4. The minimization then yields the desired distribu- 
tion g 



(5) 
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Since the problem is ill-conditioned and has many 
solutions, one single minimization does not lead to 
a reasonable (trustworthy) solution. However, we 
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FIG. 1: (a) Simulated current and the current reconstructed 
with the least-squares. The symbols (x) and (+) represent 

the real ( ) and imaginary (----) parts of the simulated 

current density j. The open symbols are the reconstructed 
currents with diflterent number of bins; (o) M/2, (>) M/4, (□) 
Af/8 and (o) A//16, where M is the number of data points, 
(b) Relative error in the logarithm of currents calculated. The 
symbols correspond to the reconstructed currents in (a) . The 
reconstructed data are calculated with M/n = 1.1 and n x 
TV ~ 10*. 



eliminate this drawback by running the minimiza- 
tion for many times with newly selected z-values. 
The resulting sets of g values are finally averaged. 
The steps 1 through 4 are executed N times with 
different sets of z (as in a Monte Carlo technique). 
The Zn and g„ values are recorded at each step, 
and the final distribution is obtained from the g{z) 
histogram. The histogram is obtained by dividing 
the 2;-axis into a number of bins (channels), and the 
average of g values is calculated in each channel, ft 
should be noted that the introduction of the Monte 
Carlo technique helps us to select independent sets 
of Zn, and the z-axis is continuous in this approach, 
which is not the case in most regularization meth- 
ods. The channel widths give the resolution of the 
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FIG. 2: Normalized polarizations g(z)/ max[go(2:)] calculated 
with the least-squares for four different numbers of bins: (a) 
M/2; (b) M/4; (c) Af/8; (d) A//16, where M is the num- 
ber of data points. The thick dashed line (----) represents 
the original polarization distribution. The distributions are 
calculated with M/n — 1.1 and n x N ~ W^. 



LIMM experiment. 

The proposed method is sensitive to the number n of 
initial unknowns and the number N of the Monte Carlo 
cycles. Best results are obtained when M/n < 1.2 and 
n X iV > 5 X 10^. It is observed that for TV > 512, the 
results do not improve significantly, and the computation 
time increases enormously. The reconstructed curves be- 
low are calculated with M/n = 1.1 and n x ~ 10*, 
which took less than 2 min on a Pentium 4 (2.4 GHz 
LinuxPc) with 1Gbyte memory. 

The original data is generated by Eq. by consid- 
ering a step polarization distribution go{z) with Az — 
4 X 10-^ 



( \ _ J rl X Az z < s/3 
^^^> ^ 1 -1 X Az z> s/3 



(6) 



We assume a 20 |J.m thick sample with thermal parame- 
ters K = 0.948 W(mK)"^ and D = 6.09 x 10"*^ m^s-^, 



3 



and A = 1 in Eq. (01 ■ The original data has 256 frequency 
points in logarithmic scale, M = 256. In the calculations, 
both the real and the imaginary parts of the current are 
used. We add both 1% Gaussian noise and a white noise 
of level of 0.1 pA to the real and the imaginary parts of 
the generated data. 

In Fig. n the simulated and the reconstructed currents 
are presented together with the relative error in the log- 
arithm of the currents. In Fig.^, the real and the imag- 
inary parts of the current density j are denoted by (x) 
and (+) symbols, respectively. The proposed method is 
very successful in restoring the original data. The data 
is reconstructed with various bin sizes. Since the imagi- 
nary part of the data has lower values than the real part, 
the relative error is larger than in the real part, and the 
restored data at high frequencies differ from the origi- 
nal data in the log-log presentation, as shown in Fig. ^p. 
Best agreement between generated and the reconstructed 
data is obtained for the largest number of bins consid- 
ered, M/2 corresponding to the smallest channel width 
2s/M. 

The normalized original and calculated distributions 
are presented in Fig. [3 In the calculations different bin 
widths are used in order to average the results. The z- 
axis is divided into 2~*M divisions with i ~ 1,...,4. 
The numerical method is able to resolve the step; how- 
ever, there are oscillations when z < s/3 which are not 
present in the original distribution. This is due to the 
applied least-squares solver, which yielded large positive 



and negative g values when two randomly selected z were 
very close. This effect could be removed by using smaller 
number of bins as presented in Fig. |21 An alternative 
way to reduce the oscilations is by implementing a con- 
strained least-squares algorithm in which the results from 
the numerical analysis with the least-squares is used as 
an a-priori assumption. In such an approach, the signs 
of the calculated g in Fig.[21are used as a pre-distribution 
for the sign of the polarization distribution in a numeri- 
cal method based on the constrained least-squares algo- 
rithm H, Q, ■ The averaging of g values is performed 
by considering various bin sizes. The obtained distribu- 
tions become closer to the original one as we increase the 
channel width (decrease the number of bins). Note that 
the resolution in LIMM is material dependent and for 
this particular generated data set, it is around 0.5 |J.ni, 
corresponding to approximately M/8 number of bins, as 
in Figs.|2t. As a results, it is not meaningful or necessary 
to use smaller channel widths to improve the polarization 
distributions. 

In conclusion, we have solved the inverse problem of 
the LIMM equation using the Monte Carlo and linear 
least-squares methods simultaneously. This method is a 
valuable alternative to those used in the literature. 

We would like to express our thanks to Dr Axel 
Mellinger for the discussions and his suggestion to gen- 
erate the data. 
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